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SOLUTION OF MATHIEUS' EQUATION ON THE ANALOG COMPUTER 



Many physical systems can be described analytically 
in terms of the functions of mathematical physics 
such as Bessel, Mathieu, and Hype rgeometric func- 
tions. However, an analytical solution of this type 
is of little value to the investigator unless it can be 
transformed into usable, numerical results. This 
transformation often is time consuming and ex- 
pensive, especially for multiple or trial-and-error 
computations. For this reason, many of these 
equations, whose practical solutions are prohibi- 
tive, are solved on a general purpose analog com- 
puter. The analog computer, which may be pro- 
grammed with ease, produces continuous, graphical 
results and allows the analyst to vary parameters 
in a few seconds for multiple computations, thereby 
reducing the time and expense required to obtain 
usable numerical results. 

Since the equations of mathematical physics de- 
scribe the behavior of a great many physical sys- 
tems, and since the analog computer is a valuable 
tool in obtaining their solution, it would be ad- 
vantageous to present the analog computer solu- 
tion to as many of these equations as possible. 
Obviously, this is impractical; however, a typical 
example can be illustrated. The illustration se- 
lected is Mathieus' equation whose solution is 
unique in that it can be stable or unstable, periodic 
or non-periodic. Mathieus' equation is a practical 
illustration, also, since it describes the behavior of 

1) wave guides 

2) moving coil loud-speakers 

3) vibrating strings and membranes 

4) frequency modulation circuits 

5) sinusoidally excited mechanical systems as 
well as other physical systems. (2) (4) (6) (7) 

This Study, then, performed on a desk-top- size 
PACE® TR-10 general purpose analog computer, 
describes the solution of Mathieus' equation. The 
objectives of the study will be threefold: first, 



Printed in U.S.A. 064 



to illustrate how Mathieus' equation should be pro- 
grammed and implemented on the analog computer; 
second, to show representative results of stable 
and unstable solutions to this equation; and third, 
to illustrate the accuracy of the computer by de- 
termining points on the stability boundary of the 
solution and comparing them to the literature 
values. 

Mathematical Model 



Mathieus' equation, which is described at length 
in several references (1, 3, 4, 5, 8), may be repre- 
sented mathematically in several forms. The form 
selected for this investigation is 



dt 



j- + (a - 2q Cos 2t) y(t) = 0 



(1) 



where y and t are dimensionless dependent and 
independent variables, respectively. The constants 
'a' 

conditions of equation (1) are 



■"' and "q" also are dimensionless. The initial 



y (o) = 1 



Idt It = 0 



o 



(2) 
(3) 



For simplicity, it has been assumed that 

a = 2q (4) 

0 < a < 5 (5) 

This restriction, which frequently occurs in practi- 
cal applications, represents the interdependence of 
system parameters upon one another and their 
practical maximum values. If one defines 



z (t) = 1 - Cos 2t 



(6) 
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then equation (1) becomes 



and 



d 2 y 

+ a z(t) y(t) =0 



(7) 



A brief discussion of the solution of Mathieus' 
equation is presented in Appendix C. 



dz 

dt z = 0 



= 0 



(10) 



This method is applicable only when the function 
to be generated is analytic, and a convenient form 
of its differential equation can be obtained. 



Computer Programming t 

As shown in Figure 1 (a mathematical block dia- 
gram of the system) the function z (t) must be gen- 
erated and interjected into the simulation to obtain 
the solution. There are two possible methods of 
generating z (t). The first method is to use a diode 
function generator which approximates the function 
over a fixed range with straight-line segments. 
This method is unsatisfactory because of the fixed 
range restriction— ±90° is typical— in addition to the 
errors introduced by the straight-line approxima- 
tion of the function. It should be noted that special 
logic circuitry can be programmed in conjunction 
with the diode function generator to provide a con- 
tinuous function; however, this only serves to point 
out the impracticality of this method. 



The maximum value of z(t) and dz/dt, which was 
determined while obtaining equation (8), is two. 
The maximum values of y(t) and dy/dt can be esti- 
mated by replacing z(t) in equation (7) by its 
maximum value to obtain a simplified version of 
Mathieus' equation, namely 



^-f + 2a y(t) = 0 (11) 
dt 

The solution to this equation— the equation of an 
oscillator— using the initial conditions defined by 
equations (2) and (3) is 

y(t) = y(o) Cos w n t (12) 



FUNCTION \ 
GENERATOR/ 



+ z(t> 



/ 



+y 



/ 



-y 



where y(o), the initial value ofy (t), is unity and the 

(13) 



frequency of oscillation, a> is 



cu = 2a 
n 



From equation (12) it is obvious that 



ay 

dt 



y(o) (Sin o) t) a> 
n n 



(14) 



azy 



Figure 1. Mathematical Block Diagram 



+2y At face value, the maximum value of y (t) is unity; 

however, since some of the solutions of interest 
in this study are unstable, y(o)— the estimated 
maximum value of y(t)— was chosen as five to 
provide a margin of safety. The maximum value of 
dy/dt then becomes twenty five, since the maximum 
value of w is 



A more accurate and efficient method of generating 
z (t) is obtained by the solution of a differential 
equation. The differential equation, which is ob- 
tained by differentiating equation (6) twice, is 



d 2 z 



y+4z(t) = 4 



(8) 



dt 



which has the initial conditions 



z (o) =0 



(9) 



max 



= ^2a 



max 



< 5 



(15) 



and the maximum value of the derivative from 
equation (14) is y (o) w^. 

Magnitude scaling is summarized in Table I for a 
± 10 volt computer. 



t It is assumed that the reader is familiar with the fundamentals of 
analog computation. 
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TABLE I - MAGNITUDE SCALING SUMMARY 



Variable 
(Dimensionless) 



Estimated 
Maximum 
Value 
(Dimensionless) 

5 

25 

2 

2 



Scale 
Factor . 
/ volts \ 
I Dim. unit / 



2 
2/5 
5 
5 



Computer 
Variable 
(volts) 



l_2y . 
'2 . 
_5 y J 
"5z ] 

[5z] 



The following scaled voltage equations* were ob- 
tained for z and y 



[5z] = 10 (J) [10] - 10 (f) [5z] 



d 
dt 



2 . 



.10 (i) ^ C 5z l 
*25 j 10 



(16) 



(17) 



The computer diagram for the simulation is shown 
in Figure 2 and the potentiometer and amplifier 
sheets, which include the static check, are shown 
in Figures 3 and 4. A tabulation of the computing 
equipment required to perform this simulation is 
contained in Appendix B. 




Figure 2. Computer Diagram 

t [ ] = reference or computer voltage 
( ) = potentiometer setting 
10, 7, etc. = input gain 



The time scale factor, /3 , was selected as one half 
so that as much of the solution as possible could 
be examined in a reasonable length of time (100 
seconds). This selection of 0 was governed also by 
the potentiometer settings which could not exceed 
one. 



PRORI PM Mathieus' Eq. 


POT 
NO. 


PARAMETER 
DESCRIPTION 


SETTING 
STATIC 
CHECK 


STATIC 
CHECK 
OUTPUT 
VOLTAGE 


SETTING 

RUN 
NUMBER 1 


NOTES 


POT 

NO. 




2/5 0 


0.800 








1 






2 


1/5 3 


0.400 








2 












3 


Constant 


0. 200 








3 


4 


1/10/3 


0.200 








4 






5 


a/5 


0.500 






Parametric Variable 


5 


6 


1/5 3 


0. 200 


0.400 






6 




7 


Constant 


0.300 








7 






a 


1/23 


0.300 


1.000 






Q 




9 


y(o)/5 


0.800 


0. 200 






9 




10 


Constant 


0.010 








10 







Figure 3. TR-10 Potentiometer Assignment Sheet 



PRORI FM Mathieus' Eq. 


AMP 
NQ 


I 

FB 


OUTPUT 
VARIABLE 


STATIC CHECK 


NOTES 


CALCULATED 


MEASURED 




OUTPUT 




OUTPUT 




' N T 


-5z 


4.00* 


-2.00 








2 


N T 


5z 


4.00 


10.00 








3 


\ 


-5z 




-10.00 








4 


\ 


-2y 




-8.00 








5 


"c 


zy 




8.00 








6 


\ 


-azy 
5 




-4.00 








7 


\ 


2 . 
" i y 


8.00 


-3.00 








8 


\ 


2y 


9.00 


8.00 








9 


\ 


t/20 


0.10 


10.00 









* 10K Feedback in Check Amplifier 

Figure 4. TR-10 Amplifier Assignment Sheet 

After examining the stability plot of the solution, 
which is derived in Appendix C and illustrated 
and tabulated in Appendix B, it was decided that 
computer runs over the range 0.5< a <^4.0 in 0.5 
increments would produce representative results. 
In addition, trial and error runs to determine the 
three transition points from stable to unstable 
solutions in this range must be made. 

Results 



Typical results of the study are shown in Figures 
5, 6, and 7. The results of all runs are summarized 
in Table II. From this summary, the stable to 
unstable transition point or stability boundary 
values of "a" for the system are 0.65, 1.75, and 
3.69 over the range of "a" investigated. These 
points are compared to literature data in Figure 8, 
which superimposes a = 2q and the computed 
values of "a" on the stability plot shown in 
Appendix B. 



TABLE II. SUMMARY OF COMPUTER RUNS 



i\umjjer 


Parameter 


Remarks 




Value 






a 




1 


0.50 


Stable 


2 


0.60 


Stable 


3 


0.62 


Stable 


4 


0.64 


Stable 


r- 

0 


0.66 


Unstable 


6 


1.00 


Unstable 


7 


1.50 


Unstable 


8 


1.74 


Unstable 


9 


1.76 


Stable 


ID 


1.78 


Stable 


11 


1.80 


Stable 


1 o 

12 


1.82 


Stable 


13 


2.00 


Stable 


14 


2.50 


Stable 


1 c 

15 


3.00 


Stable 


16 


3.50 


Stable 


17 


3.66 


Stable 


18 


3.68 


Stable 


19 


3.70 


Unstable 


20 


4.00 


Unstable 




Figure 5. y(t) versus t for a = 3.68; solution is stable on the 
verge of being periodic 
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Figure 6. y(t) versus t for a = 0.50; solution is stable, non- 
periodic 




100 120 140 160 180 200 
t 



Figure 7. y(t) versus t for a =3. 70; solution is unstable 




Figure 8. Comparison of Theoretical and Computer Results 



Conclusions 

The objectives of the study have been achieved. 
First, the mechanization of Mathieus' equation on 
the analog computer has been illustrated as well 
as several noteworthy points regarding function 
generation. A simple but accurate method of gen- 
erating an analytic function is from the solution 
of a differential equation, which generates the func- 
tion. This technique presumes that a convenient 
form of the differential equation can be obtained. 
In the case of periodic or sinusoidal functions this 
is the most practical method of obtaining a con- 
tinuous function. 

Typical solutions, which are cosine elliptic*, of 
Mathieus' equation are shown in Table I, Runs 1 
through 20 (specifically, Figures 5, 6, ' and 7). 



f Sine elliptic and cosine elliptic solutions of Mathieus 9 equation 
are defined in Appendix A. 
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These non-periodic solutions behave as expected, 
as "a" increases so does the frequency of the 
solution. By consulting the literature (5), it was 
found that the sine elliptic solutions of Mathieus' 
equation behave in a similar manner. The time 
required to obtain these results is trivial (less 
than one hour) when compared to other computa- 
tional methods. 

The percent error of the three computed stability 
boundary points compared to the literature value s(5) 
is less than 4%. This error is very small when 



one considers the error usually associated with 
the parameters used in scientific and engineering 
studies. 



The most significant result of this study is the 
fact that the stability boundary points could be de- 
termined using the analog computer. This com- 
puter application permits a system analyst or design 
engineer to select parameters and operating condi- 
tions for efficient, stable operation of a system 
with relative ease. 
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APPENDIX A 

TABULATION OF EQUIPMENT 

The following major computing components were 
required to perform this study. 

9 Operational Amplifiers 

5 Integrator Networks 

10 Potentiometers 

1 Multiplier 

1 X-Y Plotter 



APPENDIX B 



STABILITY PLOT DATA 



q 


b i 


a i 


b 2 


a 2 


b 3 


0 


1.00 


1.00 


4.00 


4.00 


9.00 


l 


-0.11 


1.86 


3.92 


4.37 


9.05 


2 


-1.39 


2.38 


3.67 


5.17 


9.14 


3 


-2. 79 


2.52 


3.28 


6.05 


9.22 


4 


-4.26 


2.32 


2.75 


6.83 


9.26 


5 


-5.79 


1.86 


2.10 


7.45 


9.24 




APPENDIX C 



ANALYTICAL CONSIDERATIONS OF MATHIEUS' EQUATION 



One of the more common representations of 
Mathieus' equation is 



2 

d y(t) 

dt 2 



+ (a - 2q Cos 2t) y(t) -0 



(1) 



where "a" and "q" are constants. The stable 
solutions to this equation, which are oscillatory, 
may be periodic or non-periodic. Fortunately, we 
need only consider those solutions which are 
periodic, since the relationships between "a" and 
"q" on an "a versus q" plot for the periodic so- 
lutions forms the stability boundaries of the solu- 
tion (5), The odd and even (sin or cos) solutions to 
equation (1) are called Mathieu functions, which 
are defined in power series as* 



ce m (t, q) = Cos m t + q c^t) + q c 2 (t) + • • • (2) 



and 



se m (t, q) = Sin m t + q s^t) + q s 2 (t) + • • • (3) 



where m denotes the order of the function. The 
" characteristic numbers" of ce m and se m are de- 
noted by a m and b m respectively (a m and b m are 
actually a in equation 1) and are related to "q" 
by a power series 

2 2 3 

a ,b = m + <* q + « q + « q + ■ • • (4) 
mm 12 3 

whose coefficients depend on the order and type of 
solution. 

An 4 'even" solution to equation (1) 



is obtained when 

y(o) = o 

and 



(9) 



m 



(10) 



A general solution to equation (1) is a linear com- 
bination of ce and se 
m m 



y(t,q) = A ce (t,q) + B se (t,q) 
' m ' m ' 

where A and B are constants of integration. 



(11) 



The coefficients of equation (4) are determined 
by substituting either equation (2) or (3) and equa- 
tion (4) into equation (1) and solving for the unknown 
c or s terms. The a coefficients of equation (4) 
are then selected to yield a periodic solution. For 
example, if m were unity 



2 



y(t).= Cos t + q c (t) + q c^t) 



dt 2 



+ q c g (t) + 



A 2 A 2 

d °1 2 d C 2 
Cos t + q — + q r-~ 



(12) 



+ q 



dt 



dt 



dt 



2 



(13) 



y(t) = ce m (t. q) 



is obtained when 

y (o) =1 

and 



(5) 



(6) 



2 

ay(t) = Cos t + q[c (t) + « 1 Cos t]+ q [c 2 (t) 

+ a 1 c 1 ( t ) + a 2 Cos t]+ • • • (14) 

and 

- (2q Cos 2t) = - q (Cos t + Cos 3t) 



fdy 

Ldt / t = 0 



= 0 



(7) 



- 2q c (t) Cos 2 t 



(15) 



while an "odd" solution 



Collecting like powers of q yields 



y(t) = se m (t,q) 



(8) 



q Cos t - Cos t = 0 



(16) 



* se and ce stand for sine elliptic and cosine elliptic. 



dt 



LIST OF SYMBOLS 



- + c ± (t) - Cos 3t + ( ot - 1) Cos t = 0 (17) 



and 



2 d2 °2 

q — — + c 2 (t) + "^(t) - 2 C]L (t) Cos 2t 



dt 



+ « 2 Cos t = 0 



(18) 



Since the particular integral corresponding to 
(a i - 1) Cos t is the non-periodic function 
1/2 (1 -a^JtSint, ^ is chosen as unity. Therefore, 



a 


Parameter 


Dimensionless 


a 

m 


= Characteristic number of ce 

m 


Dimensionless 


b 

m 


= Characteristic number of se 

m 


Dimensionless 


m 


= Order of Mathieu function 


Dimensionless 


q 


= Parameter 


Dimensionless 


y 


= Dependent variable 


Dimensionless 


z 


= Frequency variable 


Dimensionless 


A 


= Constant of integration 


Dimensionless 


B 


= Constant of integration 


Dimensionless 


& 


= Time scale factor 


Seconds 

Dimensionless Unit 


ce 

m 


= Mathieu function (cosine elliptic) 


Dimensionless 


se 

m 


= Mathieu function (sine elliptic) 


Dimensionless 



c. (t) = - - Cos 3 t 
1 o 



(19) 
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